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i . Abstract 

We present a theory of the dynamics of monatomic liquids built on 
two basic ideas: (1) The potential surface of the liquid contains three 
classes of intersecting nearly-harmonic valleys, one of which (the "ran- 
dom" class) vastly outnumbers the others and all whose members have 
■ the same depth and normal mode spectrum; and (2) the motion of 

, particles in the liquid can be decomposed into oscillations in a single 

\q | many-body valley, and nearly instantaneous inter-valley transitions 

called transits. We review the thermodynamic data which led to the 
theory, and we discuss the results of molecular dynamics (MD) simu- 
lations of sodium and Lennard-Jones argon which support the theory 
in more detail. Then we apply the theory to problems in equilibrium 
and nonequilibrium statistical mechanics, and we compare the results 
to experimental data and MD simulations. We also discuss our work 
'"O ' in comparison with the QNM and INM research programs and suggest 

directions for future research. 

O 

1 Introduction 

X 

Despite a long history of physical studies of the liquid state, no single theory 
of liquid dynamics has achieved the nearly universal acceptance of Boltz- 
mann's theory of gases or Born's theory of lattice dynamics of crystals. This 
shows the extraordinary theoretical challenge that liquids pose; they enjoy 



c3 



1 



none of the properties that make either crystals or gases relatively tractable. 
A great deal of effort has been devoted to understanding liquids as hard- 
sphere systems, which do model the core repulsion present in real liquids, 
but omit the important potential energy effects. A more realistic view was 
given by Frenkel QjJ |2] , who noted that when a typical crystal melts neither 
its specific heat, cohesive properties, nor volume changes greatly, while its 
diffusion coefficient increases dramatically; he concluded from this that the 
basic motion of particles in a liquid consists of small oscillations about a set 
of equilibria, as in a solid, but that these equilibria are neither symmetrically 
arranged in space nor unchanging in time. This highly suggestive picture 
of a liquid as something like an amorphous harmonic solid with equilibrium 
positions that occasionally move around, allowing for diffusion, has inspired 
many extensive programs of research. 

For example, after Stillinger and Weber's computer simulations [3J H] 
revealed the existence in the liquid of mechanically stable arrangements of 
particles, called inherent structures, with a wide range of energies, several 
workers have developed the idea that the liquid moves in a "rugged potential 
energy landscape" with a wide distribution of structural potential energies 
separated by barriers having a wide distribution of heights jSJIHl- These ideas 
have influenced the development of both the quenched normal mode, or QNM 
E] and instantaneous normal mode, or INM EU [TTJ H2 ECU EU ES] 
schools of thought. 

Here we will review another line of inquiry which has been pursued for the 
last five years or so [HI UH EH H3 I2DJ IHJ 122 123 121 and which differs from 
the others in that it begins by focusing on a restricted class of liquids (see 
below), and it proposes that they move in a significantly simpler potential 
landscape. 

This work is concerned exclusively with monatomic liquids, meaning ele- 
mental liquids which do not exhibit molecular bonding. Monatomic liquids 
include all elemental liquid metals and the rare gas liquids, but not molecular 
liquids such as N2 and O2, and not polyatomic systems such as alkali halides 
or water. Molecular liquids have translational, rotational, and internal vibra- 
tional degrees of freedom, while monatomic liquids have only translational 
motion, and the potential energy surface for monatomic liquids is presumably 
the simplest of all liquid potential landscapes. Our strategy is to develop a 
thorough understanding of this hopefully simplest case, and then to apply 
the insights gained there to more complex liquid systems. 

In Section |21 we describe the thermodynamic data which led us to a spe- 
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cific picture of liquid dynamics, and we describe the picture itself in some 
detail. As will be clear, the thermodynamic data are consistent with the 
picture, but they do not lead to it uniquely; thus additional support is called 
for. In Section El we review the results of molecular dynamics (MD) studies 
of two particular liquids, sodium and Lennard- Jones argon, which support 
many of our claims in far more detail. Then we apply the picture in Sections 
E] and to problems in equilibrium and nonequilibrium statistical mechan- 
ics, and we compare the results to experimental data and MD simulations. 
In Section |UJ we briefly review the picture, compare it with other research 
programs, consider the current status of its verification, discuss further prob- 
lems to which it may be applied, and describe the role we believe it fills in 
the continuing effort to develop a comprehensive theory of the dynamics of 
liquids. 



2 The picture 

2.1 Thermodynamic data 

Initial support for our picture comes from an analysis of two types of ther- 
modynamic data: The constant- volume specific heat Cy at the melting point 
of various monatomic liquids, and the entropy of melting of these elements. 

2.1.1 Specific heat 

The experimentally determined specific heats at constant pressure Cp for the 
elements have been compiled by Hultgren et al. and Chase et al. j2E] for 
both crystal and liquid phases at the melting point; these can be corrected 
in the standard way to determine Cy. Cy is composed of the contributions 
Cj from the motion of the ions and Ce from the excitation of the valence 
electrons, 

Cy = Cj + C E , (1) 

and for the nearly-free-electron elements the electronic contribution is given 
accurately by 

C E = ^ir 2 Nk 2 B Tn(e F ), (2) 

where n{ep) is the electron density of states per atom at the Fermi energy 
ep. Thus Wallace [TH] chose to study the nearly-free-electron elements, for 
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which Cj can be accurately determined. He took n{ep) from band structure 
calculations when possible [27| |2H] and from free-electron theory otherwise; 
then he subtracted out the electronic contribution to Cy, and the resulting 
ionic contributions for both the crystal and liquid phases are shown in Figure 
[TJ which is adapted from Figure 1 of [TTj- The quantities predicted by hard- 
sphere theory are shown for comparison. Notice that all elements cluster 
around Cj = 3Nks in both phases. (The exception is argon at 1 bar, which is 
known to be rather gaslike, but at pressures approaching 1 kbar its behavior 
more closely resembles that of the other liquids; thus, we will henceforth 
consider only compressed argon.) It is known that any of the crystals may be 
modeled very accurately as a set of 3N harmonic oscillators, thus accounting 
for their specific heats; this is the starting point of lattice dynamics. That 
the liquids at melt have nearly the same values for C% suggests that they too 
behave as harmonic oscillators. The departures from harmonicity for both 
phases lie outside the experimental errors; anharmonic effects on the Cj of 
liquids will be discussed in Sections EJ and El 

2.1.2 Entropy of melting 

A study of the entropy of melting at constant density (but not constant pres- 
sure) of a large number of monatomic liquids led Wallace [23 EH] to suggest 
that the elements can be separated into two classes: the "normal melting 
elements," for which AS = (0.80 ± 0.10)iV/c£, and the "anomalous melting 
elements," for all of which AS* lies far above the range of the normal melters. 
The entropy of melting results are shown in Table ^ which is adapted from 
Table III of [30J. The first column is the set of normal melters used to calcu- 
late the number given above; the second is a set of transition metals, which 
reasonably may be considered normal melters (but not anomalous melters) 
given the larger errors in their AS* data. The six anomalous melters are in the 
final column. The electronic structures of the normal melters do not change 
greatly upon melting, while for the anomalous melters the structure change 
is noticeable (semimetal crystal to metal liquid, semiconductor crystal to 
metal liquid, etc.), and this change is presumably responsible for the excess 
contribution to AS". Considering only the normal melting elements for a mo- 
ment, this change in entropy upon melting is consistent with a scenario in 
which the system, which had previously been confined to a single crystalline 
potential valley, upon melting suddenly finds itself able to move over a space 
of approximately w N valleys where Inw = 0.8; the entropy increase is due 
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C/Nk B (crystal) 



Figure 1: Ion-motional specific heat for 18 elements in both liquid and crystal 
phases at melt. Adapted from [TTj . 
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Element AS/Nk B Element AS/Nk B Element AS/Nk B 



Li 


0.75 


V 


0.90 


Sn 


1.48 


Na 


0.73 


Nb 


0.97 


Ga 


2.37 


K 


0.73 


Ta 


1.1 


Sb 


2.68 


Rb 


0.73 


Cr 


(0.9) 


Bi 


2.62 


Cs 


0.73 


Mo 


(1.2) 


Si 


3.77 


Ba 


0.90 


W 


(1.1) 


Ge 


3.85 


Fe 


0.68 


Pd 


0.74 






Al 


0.88 


Pt 


0.79 






Pb 


0.68 


Ti 


0.70 






Cu 


0.86 


Zr 


0.93 






Ag 


0.73 










Au 


0.64 










Ni 


0.88 










Mg 


0.96 










Zn 


0.97 










Cd 


0.93 










In 


0.76 










Hg 


0.90 











Table 1: Entropy of melting at constant density for 34 elements. The normal 
melting elements are in the first two columns, and the anomalous melting ele- 
ments are in the third column. Data in parentheses are less reliable. Adapted 
from [3*U] . 
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to the greater size of the available configuration space. (Strictly speaking, 
this is true only if certain other restrictions are satisfied; see Subsection 14.41 
for a more extensive discussion.) We hypothesize that this change occurs for 
all melters, normal and anomalous, and that additional factors deriving from 
the change in the electronic structure account for the difference between the 
two cases. We return to the anomalous melters in Section El 

2.2 Details of the picture 

The data considered so far suggest that particles in a liquid move in a po- 
tential landscape dominated by harmonic valleys. We have refined this ob- 
servation into a more precise picture of both the motion of the particles and 
the nature of the potential energy surface in which they move. 

2.2.1 The motion 

We hypothesize that the motion of the system may be decomposed into two 
distinct types: Oscillation in a single nearly harmonic many-body valley, and 
nearly instantaneous transitions between valleys which we call transits. That 
the valleys are nearly harmonic, and that the transits are nearly instanta- 
neous, are both suggested by the Cj data, since Ci in all cases is quite close 
to the value expected for equilibrium motion in a single harmonic valley. In 
fact, any significant departure from this behavior should show itself clearly 
in the Cj data, so we believe that this part of the picture is very solidly sup- 
ported by experiment. (Higher order corrections, to account for the small 
deviations of Cj from precisely 3Nks, are considered in Subsection 14.31 and 
discussed in more detail in Section El) We also expect transits to involve 
only a few particles in the system at a time, because transits perform a 
function in liquids analogous to that performed by collisions in a Boltzmann 
gas: They drive the system irreversibly toward equilibrium, and once there, 
they maintain equilibrium by constantly opposing fluctuations. Mechanisms 
of equilibration operate on a local level, since any small region can equili- 
brate independently of the rest (except for equilibria involving macroscopic 
coherent quantum states, not considered here), so we expect transits to op- 
erate locally. Unlike gases, though, in which collisions almost always involve 
only two particles at a time, in liquids slightly larger groups of particles can 
undergo cooperative motion, since they are sufficiently close together that in- 
terparticle potentials are always significant, so a single transit could involve 
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as many as tens of particles. 

2.2.2 Types of potential valleys 

If the valleys in the many-body potential surface are nearly harmonic, then 
each is characterized by its structure potential $ , defined to be the value of 
the system's potential energy at the bottom of the valley, and its density of 
normal mode frequencies g(u). We also hypothesize that the valleys may be 
divided into three classes: crystalline, symmetric, and random. 

(1) A crystalline valley is occupied when the system is in one of its crys- 
talline phases. These valleys are very few in number, and since the crystalline 
phases are the most stable at low temperatures, they also have the lowest 
value for their structure potential $o- Due to their very small number, the 
crystalline valleys make a negligible contribution to the statistical mechanics 
of the liquid. 

(2) The symmetric valleys correspond to more disordered configurations 
that still retain some remnant of crystalline symmetry. This group includes 
a large variety of polycrystalline and microcrystalline types, as well as the 
states of carbon realized experimentally by McKenzie, Muller, and Pailthorpe 
[3T| which differ from the perfect diamond by irregularly distorted bond 
lengths and angles. These valleys are more shallow than the crystalline ones, 
and their structure potentials $o are expected to cover a wide range of values 
due to their large variety of symmetry properties. Also because of their 
widely varying symmetries, we expect the normal mode spectrum g(u) to 
vary substantially from valley to valley. 

(3) Finally, the random valleys are occupied when the system retains no 
remnant crystalline symmetries. Since their configurations suffer no symme- 
try restrictions, these valleys should greatly outnumber both the crystalline 
and symmetric valleys; in fact, we hypothesize that almost all of the w N 
valleys available to the liquid are random, so the random valleys dominate 
the statistical mechanics of the liquid. Further, since the random valleys 
have no symmetry properties that allow them to be distinguished from one 
another, we expect that in the large- N limit all random valleys should have 
the same structure potential $ an d normal mode spectrum g(u), in stark 
contrast to both the crystalline and symmetric valleys. (In the examples we 
have studied, $o for the random valleys always lies above the <3>o values for 
all of the symmetric valleys, but we see no reason for this to be true over the 
entire potential surface.) 
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The hypothesis that the vast majority of valleys available to a monatomic 
liquid have the same depth and vibrational spectrum is a distinctive part 
of our approach, and it has extraordinary consequences for the statistical 
mechanics of the liquid; however, it is clear that the data considered to this 
point lend that idea scant support. Thus further studies were conducted to 
test the validity of this picture for specific monatomic liquids; these studies 
are discussed next. 

3 Verifying the picture 

Our picture of monatomic liquids consists of two sets of hypotheses: Those 
concerning the motion of the system, particularly that transits occur rapidly 
and involve only a few particles; and those concerning the potential energy 
surface and the classification of valleys into three types. We consider tests 
of each set of hypotheses in turn. 

3.1 Transits 

To investigate the properties of transits, we conducted computer simulations 
of two liquids: sodium and Lennard- Jones argon. 

3.1.1 Sodium 

Our simulation of an iV-atom sodium system is described in detail in [T9*] . 
The particles interact through a potential of the general form (321 122] 

mr K }) = n(V) + W <KI^ ~r L \; V), (3) 

A K,L 

where the strongly negative Q(V) is responsible for metallic binding and the 
effective ion-ion potential 0(r; V) is given by pseudopotential theory 
This pair potential is shown in Figure 1 of [19j; it is multiplied by a damping 
factor to remove long-range Friedel oscillations, and this is the only significant 
effect of the factor on the potential. After being calibrated to the bulk 
properties of crystalline sodium at K, the full potential in Equation @ has 
been shown to reproduce with remarkable accuracy several known properties 
of metallic sodium, such as the phonon frequency spectrum and the melting 
temperature as a function of pressure. In our simulations, the volume per 
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atom Va = V/N was fixed at 278 cLq, where is the Bohr radius; this is the 
density of liquid sodium at melt when the pressure is 1 bar and the melting 
temperature is 371 K. Since V is held constant in our MD calculations, we 
chose to set Q(V) to zero. The rms vibrational frequency of a typical many- 
body valley in this potential is 1.56 x 10 13 s _1 . (See Subsection 13.21 for more 
on the structure of potential valleys in sodium.) Calculations were performed 
using the Verlet algorithm [33] for a system in a cubical box with periodic 

boundary conditions; the natural time scale of the system is t* = \jlMa\l e 2 , 
where M is the atomic mass of sodium, or t* = 7.00 x 10~ 15 s. (The mean 
vibrational period in a typical potential valley is r = 57.45 t*.) The two 
parameters which varied between runs were the number of particles N and 
the MD time step St, which was always taken to be some fraction of t*. We 
will refer frequently to MD studies of sodium through the rest of the paper, 
and each time we will indicate the values of each of these parameters. 

We searched for transits in an N = 500 system where the time step was 
set to St = 0.2 1*. We cooled the system to a sufficiently low temperature that 
once it had equilibrated it remained in a single valley, as could be verified 
from its mean-squared displacement. We then raised the temperature by 
very small increments, each time allowing the system to equilibrate again, 
until transits began to occur at T = 30 K. (The details of our method of 
searching for transits may be found in [24J.) The x, y, and z coordinates of a 
particle in a typical transit are shown in Figure El Our general observations 
in sodium are as follows [21]: Every particle in the system either oscillated for 
the entire run around a single location, or it executed a transit of the general 
type seen in Figure [21 where the particle oscillated in a single region of space 
for some time, abruptly moved to a new region, and continued to oscillate in 
the new region. Typically small groups of particles transited simultaneously, 
and many more particles would execute smaller shifts in their equilibrium 
positions during a small window in time around the transit. Further, it was 
not uncommon for a single particle to participate in two or three transits, 
well separated from one another in time. The average shift in the equilibrium 
position of a particle involved in a transit was 1.75 a® (about one quarter of 
the nearest neighbor distance of 7ao); the average duration in time of any 
transit was r, and this includes the time taken by precursors and postcursors 
to some of the transits (described below in the discussion of argon). Thus our 
general picture of transits as abrupt transitions between equilibrium positions 
of a small group of particles is supported in this system. 



10 




□ 5 10 15 20 25 30 35 40 



iterations (.10') 

Figure 2: The coordinates of one particle in an 11-particle transit in sodium 
at 30 K. Adapted from |2I|. 
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3.1.2 Lennard-Jones argon 



Our simulation of argon consists of N particles interacting through a Lennard- 
Jones pair potential 



where e = 1.65 x 10~ 21 J (with equivalent temperature 119.8 K) and a = 
3.405 A. The density was set to 0.9522 particles/cx 3 , or 1.600 g/cm 3 . (The 
density of liquid argon at melt at 1 bar is 1.414 g/cm 3 .) The rms vibrational 
frequency of a typical many-body valley in this potential is 6.88 x 10 12 s _1 . 
(See Subsection 13.21 for more on the structure of potential valleys in argon.) 
Calculations were performed again using the Verlet algorithm for a box with 
periodic boundary conditions; the natural time scale of this system is t* = 



^jMa 2 /e, where M is the atomic mass of sodium, or t* = 2.16 x 10 12 s. 
(The mean vibrational period in a typical potential valley is r = 0.424 t*.) 
The time step in all MD calculations was St = 0.001 1*; the only parameter 
varied between argon runs was N. 

We searched for transits in an N = 500 system using the same technique 
that was used for sodium; we found transits at 17.1 K. The z coordinates of 
three particles involved in an 8-particle transit are shown in Figure EJ The 
horizontal dotted lines indicate the equilibrium positions of the particles be- 
fore and after the transit; the vertical line indicates the transit time. All of 
the general observations made above about transits in sodium also hold here 
pij : The type of motion seen in Figure El is typical, usually small groups of 
particles transited simultaneously, and individual particles sometimes par- 
ticipated in multiple distinct transits. The average shift in the equilibrium 
position of a particle involved in a transit was 0.44 a (about four tenths of 
the nearest neighbor distance of 1.095 a); the average duration in time of any 
transit was again r, and again this includes the time taken by precursors and 
postcursors to some of the transits. By a precursor or postcursor, we mean 
a slow drift by a single particle into a new equilibrium position either before 
or after a multiple-particle transit; a typical precursor is shown in Figure EJ 
part of the record of a 3-particle transit that occurred roughly 13 r after the 
transit shown in Figure H3 Every drift of the type seen in Figure 0] that we 
found occurred in connection with a transit, so we believe that precursors 
and postcursors are part of the transit process. They are the primary rea- 
son that the duration of a typical transit in either system is as high as r; 
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□ 2 4 6 S ID 



4 

iterations ( 10 ) 

Figure 4: The coordinates of one of three particles involved in a transit in 
Lennard- Jones argon at 17.1 K. Note the precursor, which is particularly 
visible in the x coordinate. Adapted from 
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many transits are of essentially zero duration when one neglects these effects, 
and most transits exhibit no precursors or postcursors and so are genuinely 
nearly instantaneous. Thus we conclude that our basic picture of the motion 
of particles in a liquid as a combination of oscillations and transits is veri- 
fied in these cases. The precursors and postcursors are still of some interest, 
however, and we will comment on them again in Section El 

3.2 Structure of the many-body potential landscape 
3.2.1 Sodium 

Wallace and Clements (THl ED] conducted an exhaustive study of the many- 
body potential underlying the sodium simulations in order to test the validity 
of the three-fold classification of valleys proposed above. They generated a 
large number of supercooled equilibrium states of systems with N = 500, 
1000, and 3000 and cataloged properties such as their energies and pair dis- 
tribution functions. They made the following observations about the states: 

(1) A graph of time-averaged potential energy per particle ($/N) versus 
time-averaged kinetic energy per particle (JC/N) for the equilibrium states 
is shown in Figure |SJ The melting temperature T = 371 K corresponds to 
(JC/N) = 3.53 niRy, so all of the states in the figure are supercooled, as 
claimed. This figure also shows the curve occupied by the bcc crystal states 
and the path followed by a typical MD run used to generate the states: 
Several quenches keep the kinetic energy at zero while the system moves 
down the path of steepest descent on the potential energy surface, so its 
potential energy continues to decrease; and when the quenches end the system 
equilibrates under the condition that $ + K, remains constant, so the system 
moves down the 45° line on the graph. Notice that the states separate cleanly 
into two distinct groups. Each group of states lies approximately along a 
line with unit slope, as predicted by the equipartition theorem if the states 
are moving in harmonic valleys, although the lower group shows considerable 
scatter and the slope of the upper line increases at higher temperatures. Thus 
we tentatively suggest that for each N the system is moving in a landscape 
of approximately harmonic valleys, but we also see significant anharmonic 
effects to which we will return later. 

(2) The system almost always quenched into one of the states from the 
upper group first; if the temperature was between approximately 35 K and 
200 K, it would remain in such a state for several thousand time steps (long 



15 



A 

z 

d 

V 



■10.5 



11.0 - 



11.5 



12.0 - 



12.5 



13.0 



13.5 ;i 



14.0 



■14.5 



1 1 1 


1 


/ SERIES OF QUENCHES 

/ 


m 
■ 

• 

♦ 


- 


• ♦ 


z MD RUN 


• 


/ " 


' * / 


X l/ ♦ 




\ • 

x * 


• / 






v \ :/ 








£ ♦ / 




, */ BCC STATES 


• / 


N = 500 ♦ 




N = 1000 • 


• / 

1 1 1 


N = 3000 ■ 
i 







0.5 



1.0 1.5 2.0 

< K/ N > ( mRy ) 



2.5 



Figure 5: (<&/N) versus (JC/N) for several equilibrium states in sodium. From 
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enough to compute equilibrium data) before settling spontaneously into one 
of the states in the lower group. It would remain in this state for as long as 
the MD run proceeded. 

(3) The states in the upper group lie along the same curve as the equi- 
librium states of the liquid, while the states from the lower group appear to 
be bounded in energy by the limits of the graph. 

(4) As T is increased, the graph of the pair distribution function g(r) for 
the states from the upper group smoothly evolves into g(r) for the liquid 
state, as shown in Figure |3 

(5) For a state from the lower group at temperatures above 100 K, g(r) 
exhibits a split second peak, with the first subpeak lower than the second. 

(6) By observing the mean-squared displacement of each state, defined 

by 

d(t) = -^T,[r K (t)-r K (0)}\ (5) 

Wallace and Clements found that low-temperature states from both groups 
were confined to individual valleys of the potential surface. Let d be the time 
average of d(t), or d = (d(t)) t . Then for a system in equilibrium in a single 
many-body harmonic valley, 

3h 2 T 

d= Mk^ 2 > (6) 

where G_2 (defined below) is one of the principal moments of the valley's 
frequency distribution. (For a derivation, including some subtleties involving 
the omission of zero- frequency modes corresponding to center of mass motion, 
see PH|.) Thus, if these states are confined in harmonic valleys, d should be 
a linear function of T. As we will note below, all of the valleys occupied by 
confined states in the upper group have the same frequency distribution, and 
thus the same 0_2; Figure shows d for several confined states in the upper 
group compared with Equation (JHJ) using the common value of G_2- The 
superb agreement further suggests that the valleys in which these states are 
trapped are in fact harmonic. Figure 13 of ^H] shows that the same relation 
holds for several states from the lower group, all of which are confined to the 
same valley; however, B_2 is different for different valleys occupied by lower 
states (see below), so states from different valleys fit curves with different 
slopes. 

Next, Wallace and Clements studied the actual many-body potential val- 
leys occupied by the confined states, determining properties such as each 
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2 4 6 8 10 12 14 16 18 20 

r (a ) 

Figure 6: The pair distribution function for the liquid at 390 K and pair dis- 
tribution functions for various states from the upper group at temperatures 
from 0.002 K to 201 K. Notice how the states continuously evolve into the 
liquid state with increasing temperature. Adapted from [2Hj . 



18 




5 10 15 20 25 30 35 

T (K) 

Figure 7: d versus T for several confined states in the upper group compared 
with the harmonic prediction. The inset figure shows the low-temperature 
data in more detail. Adapted from |19j . 
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valley's depth $ an d vibrational frequency spectrum g(u>). They made the 
following observations about the valleys: 

(1) The depths of the valleys occupied by the upper states all lie in a very 
narrow range, 



$ /jV = -0.01352 ± 0.00002 Ry/particle, 



(7) 



and they all have virtually the same normal mode frequency spectrum in- 
dependent of the valley or even N. The normal mode spectra for five such 
valleys are shown in Figure |H1 Since the normal mode frequencies for dif- 
ferent valleys are so similar, it comes as no surprise that the three principal 
moments of the frequency distribution, 6_ 2 , 6 , and 2 , defined by 



In k B Qo 
k B &2 



(In huj) 

~ ((M 2 



-1/2 



1/2 



(8) 



where () denotes an average over the normal mode spectrum, also vary little 
from valley to valley. (These averages always exclude the three zero-frequency 
modes that correspond to center of mass motion.) Their values fall in the 
ranges 



e_ 2 = 114±4K 
6 = 98.7 ± 0.1K 
6 2 = 154.0 ± 0.1K. 



(9) 



The larger uncertainty in 0_ 2 arises because B_ 2 is very sensitive to the 
lowest part of the frequency distribution, and thus to a small system size. 

(2) The equilibrium configuration of particles at the bottom of a valley is 
called a structure; Clements and Wallace denote the pair distribution func- 
tion for a structure G 7 (r), where 7 labels the valley in which the structure 
lies. They cooled several confined states in the upper group to find the corre- 
sponding structures and found that they all had very nearly the same G 7 (r), 
as illustrated in Figure El The fluctuations at small N gradually vanish as 
N increases. (This figure includes one valley at N = 168 which was studied 
in early exploratory calculations, but which was not used in the final work 
except at this point.) 
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200 400 600 800 1000 1200 1400 1600 



X 

Figure 8: The normal mode spectra for five different valleys occupied by 
upper equilibrium states. Eigenvalues are given as Mu\ where M is the 
atomic mass of sodium and A is a label that counts the eigenvalues. From 



21 




2 4 6 8 10 12 14 16 18 20 

r (a ) 

Figure 9: Structure pair distribution functions G 7 (r) for four different valleys 
occupied by states from the upper group. As N increases the small wiggles 
vanish. Adapted from [2*Uj . 
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(3) The universal G y (r) for the valleys occupied by confined upper states 
exhibits a split second peak (as is seen most clearly in the N = 3000 plot in 
Figure EJ), just as g(r) for the states from the lower group do [see point (6) 
above] , but with the first subpeak higher than the second. Experiments on 
Ni, Co, Cr, Fe, and Mn have identified this as a signature of an amorphous 
structure [U E3 EH1 WL ■ 

(4) Clements and Wallace also constructed the set of Voronoi polyhedra 
for each structure, and from this they computed the statistical distributions 
of two coordination numbers: The number of faces per polyhedron, and the 
angle between lines joining a particle to its Voronoi neighbors. They found 
that these distributions were universal across all of the structures found by 
cooling states in the upper group. 

(5) The valleys occupied by states in the lower group, on the other hand, 
do not exhibit universality in any of the properties measured: Their depths, 
normal mode distributions, structure pair distribution functions, and distri- 
butions of coordination numbers vary substantially from valley to valley. 

(6) The peaks in £7 7 (r) for any valley occupied by a state from the lower 
group are more numerous and narrow than the peaks in the universal G 7 (r) 
of the valleys occupied by the upper states, while the peaks in the bcc crystal 
G*bcc(^) are more narrow still. 

These results provide strong evidence that the many-body potential sur- 
face of sodium contains two distinct classes of valleys: the valleys in the first 
class exhibit universality in a wide variety of properties, while the valleys in 
the second class don't. The potential surface is dominated by valleys of the 
first class, and equilibrium states from the upper group are either confined 
to a single valley of this class or move primarily among such valleys. (That 
the first class dominates is shown by the fact that the system almost always 
equilibrates to an upper state first.) Since the upper states lie along the 
same energy curve as the liquid states and g(r) for the upper states evolves 
continuously into g(r) for the liquid as T increases, the statistical mechanics 
of the liquid is determined primarily by the properties of the first class of val- 
leys. We have less conclusive evidence that these valleys are approximately 
harmonic [points (1) and (6) about the states] and that the structures in 
the second class of valleys are less rigidly ordered than the bcc crystal struc- 
ture but more rigidly ordered than the structures in the first class of valleys 
[points (3) and (6) about the valleys]. Thus we conclude that the valleys in 
the first class are random, and those in the second class are symmetric. (The 
rms vibrational frequency and period for a "typical" many-body valley in 
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sodium in Subsection 13. II were computed for the random valleys.) Remnant 
symmetry in the lower valleys is also consistent with their lower <I>o- m this 
system, the $o values of the symmetric valleys range from the value for the 
bcc valley up to the universal value for the random valleys; symmetric valleys 
could conceivably have even higher $o values, but none were found in this 
study. Thus we conclude that our general picture of the potential surface of 
a monatomic liquid is rather well confirmed for this element. 

3.2.2 Lennard- Jones argon 

We conducted a less exhaustive study of Lennard- Jones argon [2T] in which 
we reproduced some of the results described above for sodium. One difficulty 
we found with L J argon is an interesting instability: The system has a thresh- 
old density which lies between the experimental densities of liquid Ar and 
fee crystal Ar at 1 bar, and if one attempts to cool the system at a constant 
density below this threshold, the system will collapse spontaneously until the 
threshold is reached. This limits one's ability to study densities below the 
threshold at low temperatures. Nonetheless, we were able to explore many 
valleys at the threshold density using the same technique as in sodium, and 
we found a group of states lying above those occupying the fee valley on a 
(<3>/iV) versus (JC/N) graph, as in point (1) for sodium. The valleys occupied 
by these states exhibit the following properties: 

(1) They pass both tests of harmonicity (points (1) and (6) for states of 
sodium). 

(2) Equilibrium states that move among these valleys are continuous with 
the liquid states on the (&/N) versus (JC/N) graph. 

(3) The values of &o/N for all valleys lie in the same narrow range. 
Given our experience with sodium, we conclude tentatively that we have 

found the random valleys in Lennard- Jones argon. Thus the rms vibrational 
frequency and period for a "typical" many-body valley in argon in Subsection 
13.11 were computed from these valleys. Further tests on argon and other 
liquids will be considered in Section |3 

4 Equilibrium statistical mechanics 

Now let us use the picture to develop a first order approximation to the 
statistical mechanics of a monatomic liquid. The specific heat data suggest 
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that the departures from harmonicity in the liquid's Hamiltonian may be 
treated as higher-order perturbations, at least for purposes of equilibrium 
statistical mechanics, and we shall keep this in mind as we investigate the 
Hamiltonian and compute thermodynamic quantities. Corrections beyond 
the leading order will be considered as we proceed. 

4.1 The Hamiltonian 

The general Hamiltonian for the system is written 

# = Ef§ + $({M) (io) 

Ki 

where the index K labels the particles, % labels the components of the po- 
sition or momentum of a single particle, M is the mass of one atom, and 
$ is the many-body potential. We have argued that the potential surface 
is dominated by a collection of nearly harmonic valleys; let these valleys be 
labeled with the index 7, which presumably runs from 1 to approximately 
w N . We wish to consider the form of the Hamiltonian when the system is 
localized in a particular valley. The coordinates of the particles at the valley 
bottom will be denoted {Rxi'j)}, and we define 

u K {l) =r K - R K (i) (11) 

to be the displacement of the Kth particle from its equilibrium position. The 
many-body potential in the valley will be denoted $ 7 and can be expanded 

as 

$7({«*(7)» = $0(7) + \ E ®Ki, Lj {l) «*i(7) «l;(t) + ^(7) (12) 

Ki,Lj 

where 

Ml) = = *({«Jt = 0}), (13) 

= aB^ } {{R " ]) ' (14) 

and $a(7) contains all of the higher order contributions to $ 7 . &Ki,Lj{l) is 
called the "dynamical matrix" of the potential valley. The Hamiltonian in 
the valley will be denoted H 1 and can now be written 

# 7 = $0(7) + ##(7) + M7) (15) 
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where 

H H {l) = E f§ + \ E ^, Lj (l) ««(7) ««(7) (16) 

Xi Ki,Lj 

is the harmonic contribution. An appropriate orthogonal transformation re- 
places the i*i<- (7) with new coordinates ^(7) that diagonalize the dynamical 
matrix: 

#h(7) = E (|| + q 2 x (^ . (17) 

This also defines the normal mode frequencies u}\(j). If the system contains 
N particles, then A ranges from 1 to 3iV. If the valley happens to be random, 
the Hamiltonian simplifies further; since the random valleys all have the same 
depth and normal mode spectrum, the label 7 on the frequencies and $0 can 
be dropped, so 

ff 7 = $o + ^(7)+M7) (18) 

where 

^(7) = E(^ + ^a^(7))- (19) 

The Hamiltonian in Equation (|15|) . with the harmonic part given by Equation 
(|17|) . is the starting point of our treatment of equilibrium statistical mechan- 
ics. Note that these equations describe a restriction of the full Hamiltonian, 
Equation (fTUj). to a single potential valley, so they are defined only within 
that valley. The term $^4 in the potential may describe any sort of anhar- 
monicity within the valley, but its main contribution is expected to occur at 
the edges of the valley, where the potential presumably flattens out (and de- 
parts from strict harmonic behavior) before dipping down into a neighboring 
valley. 

4.2 The partition function 

We will now compute the quantum mechanical partition function and the re- 
sulting thermodynamics, excluding exchange effects. (A quantum treatment 
is necessary for light elements including Li, Ne, and Be, but without exchange 
effects this treatment will be insufficient for describing liquid He.) We will 
also display the classical limits of our results; a fully classical development 
may be found in |IB|. 
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The partition function may be written 

Z = Tr(e-n = T, 9(E) e^ E 



(20) 



where E ranges over the eigenvalues of the Hamiltonian H and g(E) is a 
degeneracy factor which equals the dimension of the eigenspace corresponding 
to E. If the Hamiltonian described a single harmonic valley of unbounded 
spatial extent with normal mode frequencies u>\, the eigenvalues would take 
the form 

E = E (*A + ~) fiWA (21) 

where the n\ are arbitrary nonnegative integers. We have argued that the 
true potential is dominated overwhelmingly by a single class of nearly har- 
monic valleys, the random valleys, all of which have the same normal mode 
spectrum; therefore, let us approximate the eigenvalues of the harmonic part 
of the actual Hamiltonian, Equation (|17j). with values of the form given in 
Equation ()21|) using the universal random valley normal mode spectrum; 
it then remains to determine the degeneracy of each. Our approximation 
suggests the existence of eigenfunctions of H which are largely confined to 
individual valleys; clearly there would be approximately w of these for each 
eigenvalue, one per valley, and these would be approximately orthogonal 
(because they are almost spatially disjoint), hence approximately linearly 
independent. Thus we suggest g(E) « w N independent of E, or 
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(22) 



This is the approximate partition function for the liquid. In the classical 
limit {hoj\ k B T for all A), 



A 



(23) 
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We have made three noteworthy approximations in calculating Z\ First, we 
have neglected the contributions from the symmetric and crystalline valleys. 
This is a superb approximation, however, since the random valleys vastly 
outnumber the other two types. Second, we have neglected the term &a 
in the Hamiltonian (|15j). Third, in using energy eigenvalues of the form 
(J21j) we have implicitly extended a single potential valley throughout all of 
configuration space, failing to take into account its limited spatial extent; we 
have thus neglected the existence of boundaries of the valleys. These last two 
approximations are more significant, and their effects will be included in our 
subsequent calculations. 

4.3 Thermodynamic state functions 

To each thermodynamic state function X we will append a term of the form 
Xab representing the corrections due to anharmonicity and boundary effects, 
as discussed immediately above, without further comment. The Helmholtz 
free energy is 



F 



-k B T In Z 

$ - Nk B T In w + £ \-fkux - k B T \n(n x + 1)1 + F AB (24) 



where 



1 
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™A = 



e /3huj x _ i ' 



the entropy is 



S 




Nk B In w + k B l(n\ + 1) H n x + 1) - n x \nn x ] + S AB (26) 



A 



where n x is defined as before, and the internal energy is 



U 



F + TS 




(27) 
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Finally, the constant-volume specific heat is 

* ■ (§)„ 



Mn A + 1)(/3^a) 2 ] + C AB . (28) 



It is convenient to express the state functions in the classical limit in 
terms of the temperature G>o defined as in Subsection 13.21 by 

ln/c B e = — . (29) 

Using this definition, in the limit hu\ <C k B T for all A we find 

F = $o- Nk B Tlnw -3Nk B T\n(T/e ) + F AB , (30) 

S = Nk B \nw + 3Nk B [\n(T/Q ) + l} + S AB , (31) 

U = ® + 3Nk B T + U AB , (32) 

C v = 3Nk B + C AB . (33) 



4.4 Comparison with experiment 

All comparisons with experimental data will be done in the classical limit. 
First, we derive the expression for entropy of melting. The entropy for a 
monatomic harmonic crystal has the same form as the liquid, Equation (j31|) . 
without the Nk B In w term since the system resides in a single potential 
valley. Let the superscript I denote quantities of the liquid and c those of the 
crystal; then 

S l = Nk B \nw + 3Nk B [ln(T/e l ) + l} + S l AB + S E 

S c = 3Nk B [ln{T/e c ) + l]+S c A + S c E (34) 

where S B is the valence electronic contribution to the entropy (note that the 
crystal's entropy has anharmonic corrections but no boundary corrections), 
so the entropy of melting at constant density AS" is given by 

AS = S\T m )-S c {T m ) 

= Nk B lnw + 3Nk B m(e°M) + (S l AB - S C A ) + (S l E - S C E ) . (35) 
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Let us consider a normal melting element. Since its electronic structure 
is not changed significantly upon melting, it is reasonable to suspect that 
S l E ~ S E , so assuming anharmonic and boundary effects are small, the en- 
tropy of melting is dominated by the first two terms, the second of which 
depends strongly on the individual element and the first of which may or may 
not depend strongly, depending on how w varies between different substances. 
The experimental data from Subsection 12 . 1 1 21 reveal that AS* = 0.8 Nks for 
all nearly-free-electron metals with a small scatter; the only term in Equation 
()35|) that could reasonably be considered universal and thus account for these 
data is the first, assuming w is itself universal and lnu? = 0.8. That in turn 
implies that Gq ~ Q l for normal melters, with the departures from strict 
equality, along with anharmonic, boundary, and electronic entropy contri- 
butions, accounting for the scatter in AS. We have verified the prediction 
0q Oq for sodium and Lennard- Jones argon, both of which are normal 
melters; for sodium in the bcc crystal phase JH] 

6 = 99.65 K, (36) 

which is quite close to G>o = 98.7 ± 0.1 K for the liquid [Equation Q]. For 
argon, 6 = 42.5 K for the liquid at p = 1.6000 g/cm 3 , and 9 = 43.4 K for 
the fee crystal at the same density |21j . We also predict that the much higher 
AS* values of the anomalous melters can be accounted for mainly from the 
different values of 0o and Se for the two phases, both because of their very 
different electronic structures, plus the usual small anharmonic and boundary 
effects. 

Second, Wallace [TB] has compared Equations (neglecting anhar- 
monic and boundary terms) to experimental entropy data for six nearly-free- 
electron metals; the criteria used to select the six elements, and the details 
of the correction of the data for density changes, are given in ^H]. Figure 
ITU1 shows the theoretical prediction for the entropy of mercury in crystal and 
liquid phases, over a temperature range from below T m to 3.2 T m , compared 
to experimental data. The agreement is most encouraging. The differences 
between experimental and theoretical entropy as a function of T/T m for all 
six elements are shown in Figure ITT1 The differences fall within the expected 
errors in the analysis, as discussed in Recall that all thermodynamic 

functions are derived from a single potential, taken in this case to be the free 
energy, which has both a zero-temperature part, $c>; an d a thermal part (see 
Equation (J23j) or (|3*0|)): because of this, it suffices to compare one function (in 
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Figure 10: Theoretical prediction of the entropy of mercury in crystal and 
liquid phases (curve) compared with experimental data (crosses). Adapted 
from [To] . 

31 




Figure 11: Difference between experimental and theoretical entropy for six 
liquid metals. Adapted from |lfij . 
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this case, entropy) with experiment to guarantee the accuracy of the thermal 
part of our entire thermodynamic treatment. 

To check the zero-temperature part of our thermodynamics, we consider 
one further point of contact with experiment. As Wallace has shown in 
[To] , the potential minima of the liquid and crystal, corrected for density 
differences, should obey the relation 

& (Pim) ~ n(pim) « T m AS (37) 

where p\ m is the density of the liquid at melt. This assumes that both 
liquid and crystal obey the harmonic relation ($) = $o + (3/2)./V7cbT, so the 
difference in U = ($) + (3/2)NksT equals the difference in $ . However, 
as noted in ITjJ . experiments on sodium have determined that T m AS = 1.7 
mRy/atom, while MD calculations (Figure EJ) show that the difference in 
potential minima between the random valleys and the bcc valley is roughly 
0.92 mRy/atom, which is of the right order of magnitude but is 46% smaller 
than experiment. The reason for this discrepancy is easy to find; as can 
be seen in Figure and more clearly in Figure 4 of [IQ, which plots the 
states in the random valleys together with the liquid states, the slope of the 
line followed by the states increases at higher T, so (<&) does not obey the 
harmonic relation. If one uses the values of (<&) from Figure 4 of [19J at melt, 
one finds that the difference in U is in fact 1.7 mRy/atom. We will return 
to this anharmonicity in Section |U1 

5 Nonequilibrium statistical mechanics 

It is not often enough emphasized that equilibrium statistical mechanics and 
its nonequilibrium counterpart ultimately have the same starting point, the 
Hamiltonian of the system, although they make use of the Hamiltonian in 
very different ways. Both begin with a decomposition of the Hamiltonian into 
"free" and "transition" terms, where the "free" term is the more tractable 
of the two and can be fairly readily diagonalized. In equilibrium statisti- 
cal mechanics, the entire Hamiltonian contributes to the partition function, 
but often the contribution of the transition term cannot be computed ex- 
actly, so its effects are usually included as a perturbative approximation. 
Nonequilibrium statistical mechanics, on the other hand, treats the system 
as executing transitions between the states that diagonalize the free part of 
the Hamiltonian, with the transition term determining cross sections and 
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transition rates. In a gas, for example, the free part of the Hamiltonian 
describes iV-body free motion, while the transition part is responsible for 
interparticle interactions. An equilibrium treatment of a gas portrays it as 
an ideal gas with perturbations away from ideal behavior produced by the 
interaction terms; a nonequilibrium treatment via the Boltzmann equation 
portrays the gas as executing transitions between many-body free motion 
states by means of collisions, which are ultimately mediated by the same 
interaction terms responsible for perturbations in the equilibrium treatment. 
The liquid is analogous: The Hamiltonian decomposes into a "free" term, 
which is $o plus the harmonic term from Equation (jl7j) . and a "transition" 
term, which consists of $^ from Equation (|15|) plus the presence of bound- 
aries. The equilibrium statistical mechanics of the liquid, as we have seen, 
is dominated by the "free" term, although the other terms introduce cor- 
rections; and nonequilibrium statistical mechanics ultimately should treat 
the liquid as executing transits between states confined to individual valleys, 
with transits mediated by the boundary term of the Hamiltonian. (This fact 
connects transits in nonequilibrium mechanics to the boundary corrections 
in equilibrium mechanics; see Sectional) Thus transits, which do not appear 
in the equilibrium results at all to lowest order, will play a central role in un- 
derstanding the liquid's nonequilibrium behavior. Because of this we begin 
by writing the position of the Kth particle in the liquid as 

r K (t) = R K (t) + u K {t) (38) 

where -Rx(^) is the location of the center about which the particle oscil- 
lates between transits and ux{t) is motion about that center. Then Rxit) 
changes only when a transit involving particle K takes place. (Compare 
Equation This decomposition, which is motivated by a corresponding 

decomposition of the Hamiltonian, is the starting point of our nonequilibrium 
treatment. 

We will work in the linear regime, in which the coefficients determining 
the system's nonequilibrium response (self-diffusion, bulk viscosity, shear vis- 
cosity, thermal conductivity, etc.) are related to equilibrium time-correlation 
functions by expressions of the Green-Kubo form [3U|; thus our goal is to 
understand the physics behind the appropriate correlation functions. We 
will perform a sample calculation of two simple correlation functions, and 
then we will proceed to the very important velocity autocorrelation function, 
which determines the self-diffusion coefficient. We will work in the classical 
limit; quantum aspects will be discussed in Sectional 
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5.1 Correlation functions in the absence of transits 



An important part of this work involves calculating correlation functions of 
harmonically varying quantities, so we will first show such a computation by 
considering the one-particle functions (u(t) ■ u(0)) and (v(t) ■ v(0)) in the 
simplest situation, when the temperature is sufficiently low that the system 
remains in a single valley without transits. (We recall from Subsection 13.11 
that this is below roughly 30 K for sodium and 17.1 K for Lennard- Jones 
argon.) Ultimately we will be comparing these results with MD simulations, 
in which the center of mass of the system is stationary; in this case, only 
iV — 1 of the particles' positions are independent, so we define correlation 
functions as averages over particles with that restriction, and we divide by 
N—l, not N, to take into account the reduced number of independent degrees 
of freedom. We consider the position correlation function first. 

<«(*)• «(0)> = t^TT £<«*(*) -Mo)> 

iV 1 K 

= ««(°)>- (39) 

Ki 

Let the orthogonal transformation from the original coordinates to the nor- 
mal mode coordinates be denoted wja,x, so 

UKi{t) = W Ki,\ Q\(t), (40) 
A 

where the normal modes are denoted qx as in Subsection 14.11 then 
(u(t)-u(O)) = — — - £ (i»Ki,\V>Ki,\' Q\(t) qx(0)) 

Ki,X,X' 

1 



T E<»(*)gA(o)>. (4i) 

r \ 



Because the potential of the system is invariant under translations, three 
of the q\ (the components of the center of mass) are zero-frequency modes; 
since these modes are not excited by assumption, A ranges from 1 to 3iV — 3 
over the nonzero modes. Now 

(qx(t)q x (0)) = ({e- tCt qx}qx) (42) 
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where C is the Liouville operator for the system, so in our harmonic approx- 
imation 



= ( g 2 )cosKt) + ^l sinM) (43) 

IVl LOx 

where {wx} is the set of normal mode frequencies of a random valley. (The 
system is overwhelmingly likely to be in a random valley because such valleys 
dominate the potential surface.) The two averages are easily calculated in 
the canonical ensemble, 



2\ = kBT 

1x) Muo 2 ' 



) = ^T*> <<?aPa) = 0, (44) 



and the final result is 



<«W-«(0)) = — — (45) 

(Note that, aside from the fact that u(0) and u(t) would appear symmet- 
rically in the definition of the correlation function, the quantum calculation 
proceeds identically to its classical counterpart until Equation 1)440.) Since 
v(t) = u(t) in the absence of transits, a similar line of reasoning leads to 

<«(*) " «(0)> = ^^£cosM)- (46) 

A 

These results, which have the same form as the corresponding results for a 
harmonic crystal, will serve as a reference point for the work of the next 
subsection. 



5.2 Velocity autocorrelation function and diffusion co- 
efficient 

Now we will consider the velocity autocorrelation function Z(t), defined by 

Z(t) = ±(v(t).v(0)), (47) 
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which determines the self-diffusion coefficient D through the Green-Kubo 
relation [40j 

/■OO 

D = / Z(t) dt. (48) 



o 



We predict from Equation (|4l)|l that at sufficiently low temperatures 



7 rri 

Z(t) = — ~~ — %r Y cos(u x t) , (49) 
3JV-3MY 



or in terms of Z(t) = Z(t)/Z(0), 

^) = ^^E cos M). ( 5 °) 

It is not at all clear, however, how this result will be modified at higher 
temperatures by transits; certainly even an approximate solution to the sys- 
tem's equations of motion seems well out of reach. In j2Hj we argued that 
trying to understand the motion of the system in terms of a normal mode 
decomposition would be unhelpful once transits begin to occur for the follow- 
ing reasons: Over a broad range of temperatures, we suspect that any given 
particle will participate in a transit roughly once per mean vibrational pe- 
riod (this will be verified below), so each particle will experience roughly ten 
transits by its neighbors per period. Each such transit changes the many- 
body valley in which the system lies, thus changing the particular normal 
mode decomposition in which the coordinates of all the particles are ex- 
pressed. Perhaps such a change would minimally affect the coordinates of 
far away particles, but it should certainly have a substantial effect on the 
coordinates of the near neighbors. In response to this, one could instead 
suggest that the normal mode picture needs only to be supplemented, not 
replaced, and this line of reasoning has been followed most notably in some 
INM work (for example, |H1E]]); however, that work has focused not on con- 
structing an explicit model for the system's motion while transiting, but on 
modeling the effects of transits on Equation ()49)1 directly in the general form 
suggested by Zwanzig jJTj. If we must resort to models, we strongly prefer 
developing a model of the actual motion of the particles in the liquid, tran- 
sits included, and then calculating Z(t) from there, because we believe that 
the important thing to be understood is the motion, not just the behavior 
of one correlation coefficient, and because such a model can then be used 
to calculate any single-particle correlation coefficient one chooses. Thus in 
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[23*| we proposed a mean-atom-trajectory model, which consists of a single 
average particle in the liquid periodically transiting between single-particle 
equilibrium positions while executing harmonic motion between transits. We 
then incorporated into this model the essential features one expects from an 
actual solution to the equations of motion of the system, as shown below. 

Since each transit carries the system with overwhelming likelihood be- 
tween random valleys, it is sensible to model the average particle's motion 
between transits in terms of oscillations at the random valley frequency dis- 
tribution, or 

r{t) = R + u(t) 

= R + w * sin (^ + a x ) (51) 

A 

where R and u(t) are the mean- atom equivalents of Rk and Ux{t) from 
Equation (|38p. (Between transits R has no time dependence.) Now the pa- 
rameters W\ and a\ in u(t) remain to be determined. Let us assume that the 
values of the phases a\ are randomly distributed among the particles; then 
one calculates Z{t) from Equation (j51|) by differentiating to find v(t), com- 
puting the product v(t) ■ v(0), and averaging over each of the ct\ separately; 
the result is 

Z (t) = ttEI^aI^cos^). (52) 

" A 

Equation (f5*2*|) becomes Equation (f""""j) with the choice 

W " = U-lM.l W " (53) 

where W\ is an arbitrarily chosen unit vector. Thus Equation (|51|) with the 
phases a\ randomly chosen and w x given by Equation (|5*3*j). with the unit 
vectors W\ also randomly chosen, constitute our mean-atom-trajectory model 
between transits. (A brief calculation shows that this model also yields the 
correct result for (u(t) ■ u(0)) from Equation (|45p.) 

Next we must determine the effect of transits on the parameters in r(t), 
and that requires an explicit model of both the transit of an average particle 
and the rate at which transits occur. First, the transit process itself. We 
assume that the transit occurs instantaneously (the particle simply crosses 
the surface separating distinct valleys), so it must conserve both the particle's 
position r(t) and velocity v(t). To be more specific, we assume that the 
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transit occurs in the forward direction, so that the center of the new valley 
lies an equal distance away from the particle but on the opposite side from the 
center of the old valley Let r before (t), J R before , and u bclorc (t) be the position 
parameters from Equation (JoTj) before the transit, and let r after (t), i? after , and 
it after (t) be the parameters after; then our assumption of a forward transit 
implies that u ahcv (t) = -<u beforc (t), and this together with r before (t) = r after (t) 
implies 

R a{tCT = R heiore + 2u heioTe (t). (54) 

This is the change in R produced by a transit. We choose to leave the unit 
vectors W\ in Equation unaffected by transits, leaving only the effect 
on the phases a x to be determined. They must change in such a way as 
to reverse the sign of u(t) but conserve v(t); since u(t) is a sum of sines 
while v(t) is a sum of cosines, this is easily done by reversing the signs of 
the arguments {uj\t + ax) in Equation (foT|). Let the transit occur at time to] 
then uxto + af CT = -(u x t + a h x e{ore ) so 

„ after r> , . before /rr\ 

a x =-2uxt -a x . (55) 

Thus, a transit is implemented at time to by leaving the w x alone and making 
the substitutions 

R -> R + 2u(t ) 

a x -> -2Ldxt - a x . (56) 

This conserves r(t), reverses the sign of u(t), and conserves v(t). 

Let the temperature-dependent rate at which transits occur be denoted 
v, so in small time interval At a transit occurs with probability uAt. As a 
transition rate, v would ideally be calculated from matrix elements of the 
term in the Hamiltonian responsible for transits using Fermi's Golden Rule, 
and we will revisit this possibility in Section 13 but for now we will take the 
simpler path of fitting v to the results of MD simulations. 

Now the model consists of two parts, (a) Between transits, the average 
particle oscillates as given by Equations (jBT|) and with the phases ax 
and unit vectors ibx assigned randomly, (b) In each small time interval At 
a transit occurs with probability uAt; if it occurs, it replaces R and the ax 
with new values according to Equation (JHSJ). With the addition of transits, 
we can no longer express r(t) and v(t) in closed form at all times, so we no 
longer have a closed form for Z(t); but the model can be implemented easily 
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on a computer, and then the data from the run can be used to calculate Z(t) 
and Z(t) in a manner analogous to an MD simulation. 

In we calculated Z(t) in this fashion and compared the results with 
MD simulations of sodium with N = 500 and St = 0.2 1*. We performed 
equilibrium runs of the system at temperatures ranging from the glassy 
regime to nearly three times the melting temperature of 371 K. At the 
two lowest temperatures for which we ran MD, the system remained in a 
single potential valley, as could be seen from examining the mean-squared 
displacement; so these runs were compared to the model using v = 0. For 
each of the higher temperatures, we ran the model for various values of u, 
adjusting until the model matched the value of the first minimum of Z(t). 
Figures [T21 through ITB1 compare the model's predictions with a representative 
sample of our MD results; the full set of results may be found in |23j . In 
all figures, the transit rate is expressed as a multiple of r™ 1 , where r, the 
mean vibrational period in a random valley, is given in Subsection 13. II Note 
that all transit rates are on the order of r" 1 , supporting the contention made 
above that transits occur roughly once per mean vibrational period. 

The most obvious trend in Z(t) is that its first minimum is rising with 
increasing T; this is the primary reason for the increasing diffusion coefficient 
D. Note that the model is able to reproduce this most important feature 
quite satisfactorily. In fact, all fits of the model to the MD results capture 
their essential features, but we do see systematic trends in the discrepancies. 
First, note that the location of the first minimum barely changes at all in the 
model as v is raised, but in MD the first minimum moves steadily to earlier 
times as the temperature rises. The first minimum occurs at a time roughly 
equal to half of the mean vibrational period, so the steady drift backward 
suggests that the MD system is sampling a higher range of frequencies at 
higher T. Also, in Figures IT31 and IT41 the model tends to overshoot the MD 
result in the vicinity of the first two maxima after the origin, and in Figure ITBI 
this overshoot is accompanied by a positive tail that is slightly higher than 
the (still somewhat long) tail predicted by MD. These overshoots should 
clearly affect the diffusion coefficient D. To check this, we calculated the 
reduced diffusion coefficient D, the integral of Z(t), which is related to D by 
D = (ksT /M)D. The results are compared to the values of D calculated 
from the MD runs in Figure [Tol This figure includes all of the data from |2*3*] . 
including the data not reproduced here. In all of the transiting cases, the 
model overestimates D by roughly the same amount, which we take to be the 
effect of the overshoots at the first two maxima. At the higher temperatures 
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Figure 13: The model prediction for Z(t) at v = 0.35018 r 1 compared with 
the MD result for supercooled liquid Na at T = 216.3 K. From 
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Figure 14: The model prediction for Z(t) at v — 0.83985 r 1 compared with 
the MD result for liquid Na at T = 425.0 K. From [231 . 
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the discrepancy is also slightly higher, presumably due to the model's long 
tail. 

Although this single-particle model is promising, it is clearly based on 
some arbitrary choices; possible improvements, taking advantage of the in- 
formation about transits from Subsection 13.11 will be discussed briefly in the 
next section. 

6 Outlook 

6.1 What we've learned 

Two sets of experimental data on monatomic liquids, their specific heats at 
the melting point and entropy of melting, led to two hypotheses concerning 
their behavior: 

1. The many-body potential surface of a monatomic liquid is composed of 
approximately w N intersecting nearly-harmonic valleys which fall into 
three classes: crystalline, symmetric, and random. The random class 
dominates the potential surface, and in the large- N limit these valleys 
all have the same depth, vibrational spectrum, and radial and angular 
distribution functions at the valley minimum. 

2. The motion of the system decomposes into two types: Oscillation in a 
single many-body valley, and transits, which are nearly-instantaneous 
transitions between valleys. 

The picture that arises from these hypotheses has been tested successfully 
with computer simulations of sodium and Lennard- Jones argon, and it has 
been used to develop accounts of equilibrium and nonequilibrium statistical 
mechanics of monatomic liquids which compare very favorably with experi- 
ments and simulations. Both of these insights have demonstrated their value, 
and they should be taken into account in any attempt at a comprehensive 
theory of monatomic liquids. 

6.2 Further developments 

Work in any of the following areas would be of great interest. 
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6.2.1 Studies of the potential energy landscape 

The crystalline valleys have been studied for decades and are by now well 
understood. Anharmonic effects in these valleys are complicated but small 
in magnitude. We need to know more, however, about the random valleys, 
because of the dominant role they play in equilibrium statistical mechanics. 
The particle configurations of the random valley structures in sodium, for 
example, need to be characterized more completely than simply determining 
G-y(r). Do they lack any remnants of crystal symmetry, as asserted? It would 
also be worthwhile to continue the studies of argon above its critical density 
(and other noble gas liquids) until its properties are as well characterized as 
sodium's are. The remaining nearly-free-electron metals (22 or so elements) 
are expected to behave as sodium does, considering the results of pseudopo- 
tential theory for these metals; but that should be verified. Finally, most 
of the remaining elements in the periodic table that form monatomic liquids 
are non- nearly-free-electron metals (such as the transition metals), and elec- 
tronic structure theory is just now reaching the point that their interatomic 
potentials can be calculated reliably. Whether they also admit a division of 
their many-body potential valleys into similar classes would be interesting to 
discover. 

The symmetric valleys are the least studied in any of the elements we 
have considered, and they will be important if we wish to understand a real 
monatomic liquid that happens to quench into such a valley. (An experimen- 
tal example of such a case is the amorphous carbon structure in [31], cited 
in Subsection 12.21 ) The distributions of $o values and normal mode spectra 
g{oj), among other quantities, should be determined. 

Finally, we have asserted that the number of valleys is universally w N , 
where Inw = 0.8, and that the randoms so outnumber the others that vir- 
tually all the valleys are random. Is this so? Can the valleys be counted? It 
would be of tremendous interest to see if the number of valleys approximately 
obeys this relation, since it is crucial to much of the theory. 

6.2.2 Properties of anomalous melters 

Although the anomalous melters undergo substantial changes in their elec- 
tronic structure upon melting, they should obey the same liquid dynamics 
theory as the normal melters; how they become liquids should not affect how 
they behave once they are liquids. However, as we saw in Subsection 14.41 
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the values of 6 and Se for these elements should differ greatly between the 
liquid and crystal, and these differences should account for the bulk of their 
entropy of melting; testing this would be a very strong check on the theory 
of liquids we have proposed. 

6.2.3 Extensions of equilibrium theory 

We have noted that all equilibrium thermodynamic quantities have both 
anharmonic and boundary corrections, and theories of both of these need to 
be developed. Consider as an example Cj, the ionic part of the specific heat, 
which the theory predicts to be Cj = 3NIcb + Cab (cf. Equation (f3"3])). At the 
melting point, Cj for the liquid, as for the crystal, shows small anharmonic 
effects, and these appear to be of roughly the same sign and magnitude for 
both phases (see Figure |TJ). Typically the full Cy decreases as T increases, 
with Cj ultimately falling to the value for the gas, (3/2)NkB, and given 
the above comments on the anharmonic effects, we expect the boundary 
correction to be mostly responsible for this decrease. One's classical intuition 
suggests that the boundary correction is in fact negative, and this intuition 
is confirmed by calculations of the correction resulting from cutting off the 
potential of a one-dimensional harmonic valley ^7j. Further work on these 
corrections, however, remains to be done. 

Another significant anharmonic effect which is not yet understood is the 
fact that the equilibrium states occupying the random valleys in Figure 03 
and Figure 4 of ^H] do not follow a straight line of unit slope at higher 
temperatures; this is why the change in $ between crystal and liquid is not 
closer to T m AS, as discussed in Subsection 14.41 Is this feature associated 
with the onset of diffusion? Is it present in other elements, or is it unique to 
sodium? This effect is quite significant and demands further study. 

6.2.4 Extensions of nonequilibrium theory 

We have discussed only one correlation function of interest, Z(t), and it 
remains to apply the picture to several others, such as the stress-stress auto- 
correlation functions, which determine the shear and bulk viscosities, and the 
dynamic structure factor S(q, u>), which determines the liquid's neutron scat- 
tering cross section in the Born approximation. Another area of interest is 
the glass transition. It has been shown that thermal properties of a material 
during the glass transition depend on the cooling rate (321 ECU E3] and that 
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if cooling or heating stops while the system is undergoing the transition, it 
will then relax to an equilibrium state EH] • This indicates that the glass 
transition involves significant nonequilibrium behavior. A first attempt at a 
description of the glass transition using transits may be found in [22], and 
further development of that line of work is needed. 

Several questions involving the picture's conception of transits also need 
to be addressed. First, does the picture accurately portray the mechanism 
by which liquids diffuse at higher temperatures? We have seen that in low- 
temperature simulations transits occur in precisely the manner predicted 
( Subsection 13.1)1 . but that does not rule out the possibility of a qualitative 
change in behavior as temperature increases. What is the role of precur- 
sors and postcursors, which currently are not incorporated into the picture? 
Perhaps they indicate that the instantaneous transit is only a first approxi- 
mation, to be replaced by a more detailed process that unfolds over a very 
small but finite time. If the picture of transits needs to be revised, then the 
revisions should affect the nonequilibrium theory noticeably. Then there is 
the specific transit model used in our calculations of Z(t): It accounts for the 
upward shift in the first minimum in Z(t), but it requires the transit ampli- 
tude to vary as T 1//2 (because it equates the size of a transit to the amplitude 
of oscillation of a typical particle), and a softer T dependence is likely more 
accurate. Also, the transit amplitudes it predicts at the temperatures of the 
simulations in Subsection 13.11 are smaller than the observed amplitudes by 
roughly a factor of two. Then, as we have already noted, in principle one 
should be able to compute the transit rate v using the Golden Rule and the 
matrix element of the transition term in the Hamiltonian between two states 
isolated in distinct valleys. This would be a very challenging calculation, 
but it would give us tremendous insight into the mechanics of the transit 
process. Note also that v and the boundary corrections to the ther- 
modynamic quantities ultimately arise from the same source: the boundary 
term in the Hamiltonian. As such, the two should be related, and a theory 
of that relationship remains to be developed. 

6.3 The role of this theory 

This theory of monatomic liquid dynamics is based on a Hamiltonian, from 
which both equilibrium and nonequilibrium properties follow, in either quan- 
tum or classical regimes, according to the well-developed principles of many- 
body physics. The nearly harmonic character and the statistical domi- 
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nance of the random potential valleys render equilibrium statistical mechan- 
ics tractable to leading order, and they lead to well-defined corrections be- 
yond leading order. Decomposition of the motion into intra- valley oscillations 
and inter-valley transits provides a basis from which time-correlation func- 
tions can in principle be calculated from their definitions in terms of the 
mechanical motion of the system. In comparison, to our knowledge QNM 
and INM theories have been developed only for the calculation of correlation 
functions. Both work with an averaged normal mode frequency distribution 
(g(u)): QNM theories average over configurations at the bottoms of poten- 
tial valleys, while INM theories compute a temperature-dependent average 
over the entire configuration space. Although neither of these quantities 
enters the system's Hamiltonian, we can see that the QNM (g(ou)) can in 
principle approximate g(u) for a single random valley. When all is said and 
done, however, the ultimate theoretical approach to this or any other prob- 
lem is through its Hamiltonian. We believe that the ideas presented here 
provide a useful framework for thinking about monatomic liquid dynamics, 
whether one is refining equilibrium calculations to achieve improved agree- 
ment with experiment or designing and analyzing experiments to learn more 
about nonequilibrium processes. 
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